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We present a comprehensive study of the vacancy in bulk silicon in all its charge states from 2+ 
to 2 — , using a supercell approach within plane- wave density-functional theory, and systematically 
quantify the various contributions to the well-known finite size errors associated with calculating 
formation energies and stable charge state transition levels of isolated defects with periodic boundary 
conditions. Furthermore, we find that transition levels converge faster with respect to supercell size 
when only the T-point is sampled in the Brillouin zone, as opposed to a dense k-point sampling. 
This arises from the fact that defect level at the T-point quickly converges to a fixed value which 
correctly describes the bonding at the defect centre. Our calculated transition levels with 1000-atom 
supercells and T-point only sampling are in good agreement with available experimental results. We 
also demonstrate two simple and accurate approaches for calculating the valence band offsets that 
are required for computing formation energies of charged defects, one based on a potential averaging 
scheme and the other using maximally- localized Wannier functions (MLWFs) . Finally, we show that 
MLWFs provide a clear description of the nature of the electronic bonding at the defect centre that 
verifies the canonical Watkins model. 



I. INTRODUCTION 

The presence of point defects can have many conse- 
quences for the optical, electrical and mechanical prop- 
erties of materials. In particular, the behaviour of de- 
fects in semiconductors has been the subject of thorough 
and ongoing research due to their use in devices such as 
transistors, solar cells and light-emitting diodes. This 
is because point defects strongly influence the electrical 
conductivity of semiconductors by adding states in the 
band gap, thus changing the number of charge carriers 
available. Therefore, both native defects (vacancies and 
self-interstitials) and impurity-related defects (dopants) 
play a crucial role in understanding and controlling the 
performance of semiconducting materials such as silicon, 
germanium and gallium arsenide. 

To this end, the study of point defects in semiconduc- 
tors using first-principles electronic structure simulations 
has received a considerable amount of attention in recent 
yearsi Most notably, density-functional theory^ (DFT) 
has proven extremely popular due to its balance of com- 
putational speed and predictive success. Nevertheless, it 
has been noted that even for one of the simplest cases 
of point defects in semiconductors, that of the neutral 
vacancy in silicon, theoretical studies have shown a large 
scatter in results for basic quantities such as the defect 
formation energyi^ 

For a point defect in an infinite crystal lattice the re- 
laxation of atomic positions around the defect can extend 
to many successive shells of atoms. In simulations with 
period boundary conditions this gives rise to elastic inter- 
actions between the defect and its images in neighbour- 
ing supercells; the relaxation must therefore be contained 
within one supercell. Additionally, there is a spurious 
electrostatic interaction between a defect and its periodic 
replicas that depends on the size of the supercell. 

We focus on the isolated vacancy in bulk silicon (de- 
noted v 9 , where q is the charge state of the defect cen- 



tre). A number of previous studies have been under- 
taken on this defect centre using a variety of theoretical 
techniques^ in particular, Probert and Payne^, Puska et 
al&, and Wright^ have performed studies of the conver- 
gence of the defect formation energy. 

In this Article we investigate systematically and quan- 
tify the main contributions to the finite size error that 
lead to the well-known slow convergence with respect 
to system size of, for example, the vacancy formation 
energy, which arise from spurious electrostatic interac- 
tions, elastic interactions and wavefunction orthogonal- 
ity constraints between periodic images of the defect cen- 
tre in the supercell approach. In addition, our calcula- 
tions demonstrate that the defect formation energy and 
equilibrium charge state transition levels exhibit different 
convergence behaviour with respect to supercell size, de- 
pending on the Brillouin zone sampling used: the former 
benefits from the use of a dense k-point grid, while the 
latter from sampling at the T-point only. The reasons for 
this difference will be discussed in the text. Furthermore, 
we present two simple and accurate methods for calcu- 
lating the potential alignment correction to the valence 
band maximum of charged defect supercells, one based 
on averaging the potential with a Voronoi cell construc- 
tion, the other on matrix elements between maximally- 
localized Wannier functions^ (MLWFs). Finally, we re- 
late the MLWFs associated with the defect centre for 
each charge state to the canonical model of Watkins^r— 
(described in Fig. [5j Appendix [5} and show that the 
qualitative description given by this simple model is in 
full agreement with parameter-free DFT calculations. 

The rest of the paper is organized as follows: in Sec. |TT] 
we describe the computational techniques that we em- 
ploy and give the technical details of our simulations. In 
Sec. IIIII we illustrate the two methods we use to perform 
the potential alignment correction and show preliminary 
results comparing them. In Sec. IIVI we present our main 
results; first we describe the convergence properties of 
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the unrelaxed (Sec. IIV Al) and relaxed (Sec. IIVBI) neu- 
tral vacancy, and then we describe the results obtained 
for different charge states of the defect (Sec. IIV C[) in 
terms of the transition levels (Sec. IIVD[) . In Sec. IVl we 
give a brief summary of our main conclusions. 



II. COMPUTATIONAL METHODS 

In the supcrccll approximation, the formation en- 
ergy of the vacancy Ej is defined (following Zhang and 
Northrupi^) as 
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where £buik is the total energy of the bulk (perfect crys- 
tal) supercell, E q ac is the total energy of the same su- 
percell containing a single vacancy with charge q, N is 
the number of atoms in the bulk supercell, and fi e is the 
electron chemical potential. This last term can be di- 
vided into e v + A/i e , where e v is the value of the valence 
band maximum (VBM) and Afj, e is the position of the 
Fermi level in the band gap. The determination of e v is 
discussed in more detail in Sec. Mil 

For supercell sizes up to and including 256 atoms, 
the calculations arc performed using the CASTER^ 
code (version 5.0). The Ceperley- Alder local-density 
approximation^ (LDA) is used to describe exchange and 
correlation. For charged supercell calculations a com- 
pensating uniform background charge is added. Wc 
use two pseudopotentials for silicon: CASTEP's 'on-the- 
fly' Vanderbilt ultrasoft pseudopotentiat^, and a norm- 
conserving pseudopotentiaU% both have four valence 
electrons and give accurate lattice constants for bulk sil- 
icon (5.39 A and 5.38 A respectively, compared with an 
experimental value of 5.43 A). This slight underestima- 
tion of less than 1% is typical for LDA DFT. The DFT- 
optimized lattice parameter is used in all calculations. 

For our largest calculations, on a 1000-atom simple 
cubic (SC) supercell, we use the linear-scaling DFT code 
ONETE P 17 ' 18 (version 2.4). The 'cross-over', namely, the 
system-size at which it becomes computationally more 
efficient to use ONETEP as opposed to conventional cubic- 
scaling DFT, is highly system dependent and lies at 
around 500 atoms for silicon. For this reason we use 
conventional plane- wave DFT for supercells smaller than 
this cross-over, and linear-scaling DFT for the largest su- 
pcrccll. ONETEP makes use of the single-particle density 
matrix, expressed in separable form in terms of a local- 
ized basis of non-orthogonal generalized Wannier func- 
tions (NGWFs) and a density kernelJi^ We use the 
same norm-conserving pseudopotential as for CASTEP, 
and nine NGWFs on each silicon atom with a trunca- 
tion radius of 3.97 A. We do not truncate the density 
kernel. 

We follow the methodology outlined by Probcrt and 
Payne in their study of the neutral silicon vacancy^; how- 
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FIG. 1. (Color online) Potential alignment correction for v 2_ 
using the MLWF and real-space Voronoi cell methods. The 
values are adjusted with respect to the bulk supercell. The 
dashed horizontal lines show the value for the furthest point 
from the defect centre. 



ever, numerical parameters are converged independently 
for all system sizes. The supercells are constructed from 
three unit cell shapes: FCC (the primitive cell, with 2 
atoms), SC (with 8 atoms), and BCC (with 32 atoms). 
The supercells are then made from n 3 unit cells. We 
perform calculations on both unrelaxed and relaxed ge- 
ometries; our convergence tolerance for the defect for- 
mation energy of a given supercell is 10 meV. We use 
a regular Monkhorst-Pack (MP) mcshi^ of k-points for 
the Brillouin zone integration. In the text, the param- 
eter /cmp refers to a k^p grid. The issue of Brillouin 
zone sampling shall be discussed in Sec. IIVI in general, 
we make use of two sampling schemes: a T-point only 
sampling (i.e., &mp = 1); an d what we call a 'dense' 
sampling. This last term we shall use to indicate that 
&mp has been converged with respect to the formation 
energy for a particular supercell. 

For the relaxation a quasi-Newton BFGS scheme is 
used. All the atoms in the supercell are allowed to move, 
and their positions are slightly randomized at the start 
of the procedure to allow for symmetry breaking. Our 
convergence tolerance for the geometry optimization is 
5 x 10~ 3 cV/A for the root mean square force on all the 
atoms. 



III. DETERMINING THE ELECTRON 
CHEMICAL POTENTIAL 

As described in Sec. [HI the formation energy of charged 
defects depends on the electron chemical potential, which 
is given relative to the VBM eigenvalue position e v . e v is 
more precisely defined as the energy difference between 
the pure host material and the host with a single elec- 
tron hole. In the limit of a dilute hole gas, this energy 
difference is equivalent to the VBM eigenvalue. For a fi- 
nite supercell calculation, however, the energy difference 
between the bulk system with and without the hole only 
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slowly converges to the infinite limit as the system size 
is increased^ Instead, the VBM will always remain at 
the same position (given a sufficient Brillouin zone sam- 
pling), which is the limiting value. 

Two approaches can be taken to determine e v : either 
calculating the energy difference of the bulk supercell 
with q = and q = +1, or using the bulk VBM eigen- 
value and employ a potential alignment correction^— . 
In the rest of this section we describe the two methods 
that we employ for calculating the potential alignment 
correction: a partitioning of the real-space potential us- 
ing Voronoi cells, and an MLWF-based approach. To the 
best of our knowledge, this is the first time these methods 
have been used to calculate the correction. 



A. The real-space Voronoi cell method 

Our first method is based on considering the electro- 
static potential directly. In such approaches, the correc- 
tion is typically determined by plotting the electrostatic 
potential obtained from the DFT calculation, and aver- 
aging it either in the x-y plane, within atomic spheres 
or over primitive cells of the host material. The impor- 
tant point is that the averaging must be done in some 
localized manner, such that it is possible to measure the 
value of the bulk-like potential in the charged defect su- 
percell by considering a region far away from the defect 
centre. The correction to the defect formation energy is 
then q (V va _ c — Vtmik), where Vt> u ik and V vac are the aver- 
age potentials in the bulk supercell and in a 'bulk-like' 
region of the vacancy supercell. 

In order for the correction to be unambiguous and ac- 
curate, the averaging volume should be as small as pos- 
sible while still covering a region of space that is com- 
pletely representative of the bulk material (i.e., all re- 
gions in the bulk supercell should give the macroscopic 
average for the material) . Our approach is to use Voronoi 
cells: each region corresponds to the real-space volume 
of points closest to one particular atom in the supercell. 
There are, therefore, as many regions as there are atoms. 
The electrostatic potential, however, is given on a dis- 
crete real-space grid whose spacing depends on the basis 
set cut-off energy and which, in general, is not commen- 
surate with the atomic spacing, leading to inaccuracies in 
the averaging procedure that may swamp the difference 
between the bulk and defect supercells. In order to over- 
come this drawback, we Fourier-interpolate the potential 
onto a finer grid that is commensurate with the atomic 
spacing, which results in each Voronoi cell containing ex- 
actly the same number of points. We note that care must 
be taken to give the correct fractional weighting to points 
which arc directly between two or more atoms. 

The potential alignment correction is then determined 
by considering Voronoi cells belonging to the atoms which 
are furthest from the defect centre, as shown in Fig. [TJ 



B. The MLWF method 

Our second approach to calculating the potential align- 
ment correction is based on matrix elements of the Hamil- 
tonian in the basis of MLWFs. The shift in the reference 
of potential that is sought is reflected in the position of 
the energy eigenvalues; however, comparing the eigenval- 
ues of bulk-like states in the defect supercell with those 
of the bulk supercell is problematic since the eigenstates 
are delocalized. MLWFs, on the other hand, provide a 
probe of the localized properties of the system. Once 
the Hamiltonian is transformed into a basis of MLWFs, 
the potential alignment correction can be determined by 
considering on-site matrix elements H nn = (uj n \H\uj n ) 
(where u> n is a MLWF belonging to the cell at the ori- 
gin) for Wannier functions whose centres, defined by their 
first moments, are far away from the defect centre. The 
shift in potential is then simply H nn ^ ac — -ff mi! buik- In 
practice, as shown in Fig. [I] we plot this difference for all 
n as a function of the distance of the Wannier function 
centre from the defect centre, and the value of the po- 
tential alignment correction is taken as the mean of the 
values which are furthest from the defect centre. 

This method is feasible since the wannierization pro- 
cedure results in the same partitioning of the electronic 
density in the bulk and defect supercells everywhere ex- 
cept in the direct neighbourhood of the defect; hence it 
is possible to identify Wannier functions that are associ- 
ated with the defect centre and Wannier functions that 
are not. In the latter case, therefore, it is straightforward 
to unambiguously match equivalent Wannier functions in 
the bulk and defect supercells. 

From Fig. Q] it can be seen that the MLWF correc- 
tion method is in excellent agreement with the Voronoi 
cell method. Among all our calculations for the various 
charge states, the maximum discrepancy was 0.01 eV. 
The MLWF method gives a finer representation of the 
system than averaging over atomic sites (as in the 
Voronoi cell approach), since each silicon atom has four 
bonding Wannier functions connecting it to its neigh- 
bours. Therefore, there are twice as many Wannier func- 
tions as atoms, and hence twice as much information to 
consider^ 

MLWFs are now a standard tool for the analysis of 
electronic structure calculations; therefore, this approach 
provides a simple and accurate method for calculating 
the potential alignment correction, since all the neces- 
sary information is readily available once the standard 
wannierization procedure has been performed^ 

IV. RESULTS 

A. The neutral unrelaxed vacancy 

In order to study the finite size convergence proper- 
ties of the system we first simulate the neutral vacancy 
in its unrelaxed state, with all the atoms in their perfect 
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TABLE I. List of all supercells up to 256 atoms with their 
respective symmetries. Also listed are the value of /cmp used 
(converged with respect to the formation energy for each su- 
percell), the corresponding k-point volume in reciprocal space, 
the unrelaxed defect formation energy and the kinetic energy 
contribution to this value. Calculations were performed with 
an ultrasoft pseudopotential and a plane-wave energy cut-off 
of 400 eV. 

crystalline positions; our results are reported in Table HI 
Considering the different supercell geometries separately, 
it can be seen that increases monotonically with sys- 
tem size. We therefore take our highest value of 4.13 eV 
(for the 256-atom BCC supercell) as a lower bound to 
the unrelaxed defect formation energy. 

The slow convergence of the formation energy as a 
function of system size is due to the spurious interac- 
tion between periodic images of the defect centre. As the 
defect is both unrelaxed and uncharged, this is neither 
caused by elastic interactions nor monopole-monopolc 
electrostatic interactions; rather, it is the result of higher 
multipole interactions and overlap between the wavefunc- 
tions of the periodically repeated defect centres. 

The importance of wavefunction overlap may be 
gauged by considering the kinetic energy contribution to 
the defect formation energy, which we define analogously 
to Eq. <[T|), except using only the non- interacting kinetic 
energy contribution to the total energy instead of the to- 
tal energy itself. The results are listed in the right-hand 
column of Table HI It is immediately apparent that the 
kinetic energy component of the defect formation energy 
varies on a scale larger than the total formation energy 
as the supercell size is increased, and that, even at our 
largest supercell sizes, it is not converged to better than 
0.1 eV. These changes are caused by subtle changes in 
the defect wavefunctions of neighbouring defect centres 
that occur in order to maintain orthogonality as the pe- 
riodic images of the defect move apart. While the elec- 
tronic density is rather insensitive to them and, therefore, 
the electrostatic interactions are not affected, the kinetic 
energy is not. As a result, any correction scheme that 
accounts solely for the classical electrostatic interaction 
will not be sufficient in this case to predict the correct 
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FIG. 2. (Color online) Dispersion of the unrelaxed defect 
level. The symbols show the position of the level at T, and 
the bars show the extent of the dispersion obtained from our 
dense k-point sampling. The positions are given with respect 
to the bulk VBM. 

formation energy for an infinite system beyond this level 
of accuracy. 

We consider one last question before moving on to 
ionic relaxations: is it desirable to use a dense k-point 
grid when simulating point defects? Although this is 
obviously necessary for describing the delocalized bulk 
states, it seems reasonable that a T-point sampling might 
give a better approximation of the localized defect levels 
in the dilute limit, as a way of deliberately eliminating 
unwanted dispersion, and by occupying the state which 
most closely resembles the bonding defect state for the 
infinite system^. Indeed, Fig. [5] shows that the position 
of the defect level at T rapidly converges to a fixed value, 
even if the total dispersion of the level throughout the 
Brillouin zone is large. 

As far as the defect formation energy is concerned, at 
least, our results show that T-point calculations do not 
give a better estimate of the defect formation energy: for 
the 256-atom supercell E® is underestimated (as com- 
pared to our lower bound estimate) both for the unre- 
laxed and relaxed cases, by 0.2 eV and 0.4 eV respec- 
tively. This conclusion is in agreement with two other 
studies on the relaxed silicon vacancy: Puska et al&, who 
report that a 2 3 MP sampling shows a faster convergence 
than a T-point sampling for the neutral charge state; and 
Wright^, whose results for the 216-atom supercell show a 
greater average error across all charge states for T-point 
sampling than any finer MP grid sampling. We note 
that, conversely, quantities of interest other than the de- 
fect formation energy might benefit from a T-point only 
sampling; as we shall see in Sec. IIVD1 the stable charge 
state transition levels are an example of this. 



B. The neutral relaxed vacancy 

On relaxation of the ionic positions, the predicted sym- 
metry for v° from Watkins' model (D2d) is not seen in 
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FIG. 3. (Color online) Relaxation effects for the neutral va- 
cancy. Squares show the defect formation energy for the un- 
relaxed lattice and circles for the relaxed one. Labelled below 
is the point group symmetry of the defect centre after relax- 
ation. The last column shows the relaxed formation energy 
from previous studies: a) Ref. 0] (N = 216, fc M p = 2, LDA); 
b) Ref. @| (N = 216, fc M p = 2, GGA); c) Ref. [| (N = 256, 
k M p = 2, GGA); d) Ref. {N = 1000, fc M p = 3, LDA); e) 
Ref. 0| (N = 1000, fc M p = 3, GGA). Our calculations were 
performed with an ultrasoft (norm-conserving) pseudopoten- 
tial and the LDA functional, with a converged plane-wave en- 
ergy cut-off of 400 eV (800 eV) and a k-point grid of &mp = 3 
(k M p = 4). 



any supercell smaller than 256 atoms, as shown in Fig. [3] 
Calculations on the 32-atom BCC, 64-atom SC and 128- 
atom FCC supercells do not exhibit a change in symme- 
try, only an inwards relaxation. Consequently, the defect 
formation energy is only lowered by ^0.2 eV in these 
cases. 

The 256-atom BCC supercell undergoes the predicted 
change in symmetry and a reduction in defect volume 
by more than 40%. This results in the defect formation 
energy being lowered by 0.7 eV (full results are given in 
Table ID). Therefore, even though the unrelaxed defect 
formation energy is highest for the largest supercell, af- 
ter relaxation it becomes the lowest, demonstrating the 
well-known importance of Jahn- Teller distortion and the 
long-ranged nature of the clastic interactions which can 
lead to both qualitatively incorrect relaxation patterns 
and quantitatively inaccurate formation energies in small 
supercells. 

These results confirm several previous DFT studies on 
this systemi££ In particular, Puska et alfl report that 
small supercells can show a range of different symmetries 
depending on the k-point sampling, and there is some 
evidence from their results that a T-point only sampling 
favours D2d- It is clear from our results that a dense k- 
point sampling favours the lattice's unrelaxed tetrahcdral 
point group symmetry Tj for small supercells. 

Fig.[3]also shows the relaxed defect formation energy in 
a 1000-atom SC supercell, as calculated with the linear- 



scaling DFT code ONETEP<2& The result (calculated with 
the same pseudopotential) is within 23 meV of the 256- 
atom BCC system, and the relaxed defect structure is 
almost identical, with only a 33 mA root mean square 
difference in the lengths of the bonds at the defect centre. 



C. Charged vacancies 

Our results for the relaxed symmetries and defect for- 
mation energies of the various charge states of the va- 
cancy are shown in Table [TTJ All the charge states apart 
from v 2 ~ are in agreement with Watkins' model. For 
v 2 ~ a completely different 'split vacancy' configuration is 
found, with a point group symmetry of D3d (in agreement 
with previous studies by Nieminen et al£&- and Wright^) . 
The symmetry predicted by the model (C2 V ) was found 
as a metastable state. For all charge states there were 
no qualitative differences between the T-point only and 
multiple k-point calculations, although the former con- 
sistently showed a greater inwards relaxation and thus a 
smaller defect volume. For a detailed discussion of the 
electronic structure of the different charge states of the 
vacancy, the reader is referred to Appendix IA1 

We note that the relaxation of v 1_ is the most prob- 
lematic; in fact, if the initial random displacement of the 
atomic positions is not large enough, the system consis- 
tently converges to a metastable Td configuration with 
a completely spin-polarized triply degenerate defect level 
(i.e., filled with three spin-aligned electrons). Further- 
more, our final lowest configuration is only approximately 
C2v, as there is a small but noticeable splitting of ~0.1 A 
of the bond lengths between the two pairs of neighbours 
of the vacancy. Nieminen and Wright both report that 
the LDA does not give the expected symmetry for this 
charge state, although they find a D3d configuration. 

As explained in Sec. IIII1 there arc two possible ap- 
proaches for estimating the position of e v in the system: 
we can either use a calculation of the bulk supercell with 
an electron hole (denoted 'Hole' in Tables HTMlIIf) . or the 
bulk VBM value with a potential alignment correction 2 ^ 
('VBM'). The agreement between the two approaches is 
better for the positively charged vacancies (—0.05 eV to 
0.01 eV) than the negatively charged ones (0.05 eV to 
0.21 eV); however, these uncertainties are small enough 
not to affect the level ordering. 



D. Transition levels 



Following the careful description of Baraff et al*^, we 
define the stable charge state transition level E (m/n) to 
be the value of the Fermi level (with respect to the perfect 
crystal valence band edge) at which there is a crossing of 
the defect formation energies of two charge states m and 
n, leading to a change in the most stable state from v m 
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For this system we use the relaxed atomic positions obtained for the 1000-atom v° supercell. without further relaxation. 



Bond 



TABLE II. Summary table of the results obtained for the 256-atom BCC supe rcell a nd the 1000-atom SC supercell 
lengths are given between the atoms surrounding the vacancy, as labelled in Fig. 6(a) for v 2+ , v 1+ , v°, V 1 " , and Fig. |6(g)| for 
v 2 ~. Atoms labelled with the same letter are equivalent by symmetry. Aa = a — 90 is the distortion of the regular octahedron 
in the split vacancy configuration. The defect volume is calculated from the tetrahedron formed by the four neighbours of the 
vacancy site. The defect formation energy is given for Afj, e = 0. Symmetry labels are given in Schonflies notation. 

















Ref. [6] 


Ref. [3£ 


Ref. [7] a 
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N = 


256 




N = 1000 


N = 216 


N = 216 


N = 


1000 




fcMP 


= 1 


&MP 


= 3 


k-MP — 1 


kup = 1 


&MP = 2 


&MP 


= 3 






Hole 


VBM 


Hole 


VBM 


VBM 


LDA 


GGA 


LDA 


GGA 


E(2+/l+) 


0.13 


0.13 


0.16 


X 


X 




0.19 


0.13 


0.27* 


0.27 


E{2+/0) 


0.09* 


0.10* 


0.11* 


X 


X 




0.15* 


0.10* 


0.28 


0.19* 


E(l+/0) 


0.05 


0.07 


0.06 


X 


X 


0.04 


0.11 


0.06 


0.28* 


0.12 


£(o/i-) 




0.43 


0.38 


0.35 


0.26 




0.57 


0.37 


0.76 


0.63 


£(0/2-) 




0.32* 


0.25* 


0.24* 


0.13* 




0.49* 


0.37* 


0.66* 


0.53* 


£(W2-) 




0.22 


0.13 


0.13 


0.00 




0.40 


0.36 


0.56 


0.42 



a The transition levels have been calculated from the quoted values of the formation energy of the various charge states using Eq. J2} . 



TABLE III. Transition levels for the 256-atom BCC supercell and the 1000-atom SC supercell. Results from previous studies 
are also shown. Our calculations use the LDA functional. Crosses (x) are used when there is no transition in the band gap. 
Asterisks (*) denote thermodynamically stable transitions. All values are given in eV. 



to v" as the Fermi level is raised. This is given by 

E n f (Au e =0)-Ef (Au e = 0) 

E(m/n)= f } -. (2) 

m — n 

As first predicted by Baraff et al^- and then experi- 
mentally verified by Watkins and Troxell 2 ^, the isolated 
silicon vacancy exhibits a 'ncgative-U' effect, by which 
the stable charge state changes directly from v 2+ to v°. 
The v 1+ state is therefore only metastable; this is a direct 



consequence of the Jahn- Teller distortion, which lowers 
the energy of the neutral vacancy much more than the 
singly positive vacancy. It has also been suggested that 
the same effect might be observed in the sequence v°, 

Our transition levels are given in Table IIIII as can be 
seen, the T-point calculations predict the negative-U be- 
haviour both for the negatively and positively charged 
levels as expected, and are in good agreement with the 
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Stable charge state 
Expt. levels 







1 


1 1 1 1 














1 


i i i i 



0.1 0.2 0.3 0.4 0.5 

(eV) 



FIG. 4. (Color online) Formation energy of the different 
charge states of the vacancy as a function of the electron 
chemical potential (plotted relative to the VBM). The en- 
ergy range shown covers the DFT band gap for silicon. The 
thermodynamically stable charge state at each point is high- 
lighted in bold, and the circles indicate the level position for 
the stable transitions. The dashed vertical lines show the 
available experimental values for the transition levels. The 
results shown are for the 256-atom BCC supercell with a F- 
point sampling. 

available experimental results (also shown in Fig. 2]). Sur- 
prisingly, the calculations using a dense k-point sampling 
= 3) give the opposite ordering for the sequence 
v 2+ , v 1+ , v°; this results in no transitions in the band gap 
between these levels, and only v° being a stable charge 
state. The negative-U behaviour is however still present 
for the negatively charged levels, although their position- 
ing is lower than for the T-point calculations. The T- 
point, therefore, gives better estimates of the transition 
levels. However, as stated previously, the absolute value 
of the defect formation energy is not as well converged 
with respect to the dilute limit when using this Brillouin 
zone sampling. However, calculations of v° and v 1+ in 
the 1000-atom supercell at T (also included in Tables |HT - 
IIIip are converged with respect to system size both in 
terms of the transition level E (1+/0) = 0.04 eV and the 
absolute value of the formation energy E® = 3.55 eV. 
Therefore, this suggests that at this system size a T-point 
sampling can be employed for simultaneous convergence 
of the both quantities of interest. 

V. CONCLUSIONS 

We have studied the silicon vacancy in all its charge 
states using the supercell approach and plane- wave pseu- 
dopotential DFT. Our calculations confirm the slow finite 
size convergence of defect formation energies and transi- 
tion levels, due to electrostatic interactions and wave- 
function overlap between periodic images of the defect, 
and long-ranged ionic relaxations. The impact of each of 
these has been quantified, and it has been found that all 
three provide non-negligible contributions to the total er- 



ror. Furthermore, due to the hybridization of the defect 
levels with the perfect crystal band structure, we find 
that the choice of k-points has a noticeable impact on 
the results. In particular, T-point calculations converge 
faster than calculations with uniform, multiple k-point 
sampling when considering stable charge state transition 
levels (given by differences in formation energies of differ- 
ent charge states), and vice- versa when considering the 
absolute value of formation energies. 

Due to this slow convergence, the calculations may still 
benefit from the use of even larger supercclls than the 
ones currently employed. This, however, will present ad- 
ditional numerical challenges: the defect formation en- 
ergy is an intensive quantity that is obtained by taking 
the difference of two total energies that are extensive with 
the system size. As a result, to achieve a given accuracy 
for Ef as the supercell size is increased, it is necessary to 
proportionately increase the precision per atom in the to- 
tal energy calculations in order to avoid the result being 
swamped by numerical noise. 

We have also introduced two methods for correcting 
the alignment of the valence band maximum for charged 
defects: one is based on averaging the electrostatic po- 
tential using a Voronoi cell construction, and the other 
on Hamiltonian matrix elements in a basis of maximally- 
localized Wannier functions. The two methods give ex- 
cellent mutual agreement and constitute simple and ro- 
bust ways to calculate potential alignment corrections. 
However, the determination of e v remains a somewhat 
uncontrolled source of error, as demonstrated by the 
discrepancy between the potential alignment correction 
('VBM') and total energy difference ('Hole') approaches 
(described in Sec. IIIip . which can be as small as 0.01 eV 
(for positively charged defects), or as large as 0.2 eV 
(for negatively charged defects), as shown in Table |TT] 
(right hand columns). We note that the problem of cor- 
recting errors in charged supercclls arising from periodic 
boundary conditions has been addressed by many previ- 
ous studie a 21 ' 29 " — . 

The accuracy of the LDA functional for describing ex- 
change and correlation should also be investigated. GGA 
calculations have already been shown to give very similar 
qualitative and quantitative results (as shown in Sec. HVP ; 
both local and semi-local functionals, however, suffer in 
particular with respect to the well-known effect of self- 
interaction error. Unfortunately, higher accuracy meth- 
ods are at present confined to relatively small system 
sizes due their computational expense, and are therefore 
limited by the large finite size errors that we have dis- 
cussed earlier. Nonetheless, qualitative differences have 
been suggested for the vacancy from screened-exchangc 
calculations^ 

Finally, it is interesting to note that maximally- 
localized Wannier functions give a chemically intuitive 
picture of the electronic structure at the vacancy site, 
and without any external input confirm Watkins' LCAO 
model of the defect. The only charge state which does not 
follow the predictions from the model (the doubly neg- 



ative vacancy) is also explained in terms of its Wannicr 
functions. Full details are given in Appendix |A"1 
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Appendix A: Visualization using MLWFs 

Wannicr functions provide a localized view of the elec- 
tronic structure. This is particularly useful in the study 
of point defects, since the Bloch states associated with 
the defect levels are not completely localized due to the 
periodicity of the system; in fact, in small supercells it 
is not possible to disentangle them from the bulk band 
structure. 

In the case of bulk silicon, previous work has shown 
that it is possible to recover the typical a bond orbitals 
between silicon ions by wannierization of the valence 
band; alternatively, the bottom conduction bands (once 
disentangled from higher bands) give the corresponding 
antibonding orbitals, while treating both these sets of 
bands together as a single manifold produces four sp 3 or- 
bitals on each ion£i2£ The bulk silicon a orbitals obtained 
from our calculations are shown in Fig. |6(b)[ 

For the wannierization of the defect systems we only 
use the occupied manifold, except in the case of the unre- 
laxed vacancy where we include all the defect levels. We 
can loosely define the concept of a defect Wannier func- 
tion to be any MLWF in the defect supercell which differs 
qualitatively from the a bond MLWFs of the bulk sys- 
tem. For all charge states, such defect Wannier functions 
are only present within the first ionic shell around the 
vacancy; between the first and second shell we already 
recover cr-like bonding orbitals. 

Fig. 6(c) shows the defect Wannier functions for the 
neutral unrclaxcd system; as expected from Watkins' 



model (described in Fig. [5]), these are sp 3 orbitals point- 
ing towards the vacancy. These orbitals cannot be ob- 
tained by wannierizing the three visible defect levels in 
the gap alone, as the nodeless combination which lies 
within the valence band is also needed; therefore, the 
entire valence manifold plus the defect levels in the gap 
must be wannicrized as a whole. 



v 



2 + 



V 



II 



The defect Wannier functions obtained for 

in their relaxed configurations (Figs. 6(d) 6(f)) 



and v 

also support Watkins' model. For the doubly positive va 
cancy, the defect levels in the gap are empty and there 
fore not included; the result is a single MLWF corre 
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FIG. 5. (Color online) Watkins' LCAO model for the 
silicon vacancy, deduced from electron paramagnetic res- 
onance (EPR) studies of the defect centre—, and later 
confirmed by electron-nuclear double resonance (ENDOR) 
measurement :] 37 ' 38 . The four orbitals associated with the 
neighbouring atoms of the defect site are the dangling bonds 
resulting from the removal of the central silicon atom. The fig- 
ure shows the predicted ionic configuration and point group 
symmetry for different charge states of the vacancy due to 
Jahn- Teller distortion 3 ^. For the D2d and C2v configurations, 
the first charge state (in black) refers to the spin-polarized 
case with only one electron in the highest occupied orbital, 
and the second charge state (in red) refers to the non-spin- 
polarized case with two electrons. 



sponding to the symmetric s-type nodeless combination 
of the four sp 3 orbitals. The neutral vacancy includes one 
lowered defect level only, and produces the two bonding 
orbitals between the pairs of nearest neighbours of the 
defect centre. The singly negative vacancy includes two 
defect levels in the gap (for the majority spin), produc- 
ing one bonding orbital (between the pair of ions with 
the shortest bond length) and two sp 3 orbitals. 

Finally, we show the split vacancy configuration ob- 

In this case, one of 



tained for v 2 in Figs. 6(g) 6(h) 



the neighbours (labelled c) moves halfway along the line 
connecting it with the vacancy site, thus placing itself 
at the centre of an octahedron made up of two pairs of 
second- nearest neighbours (labelled a and b). This octa- 
hedron is approximately regular (having as faces equilat- 
eral triangles); however, there is a small distortion due 
to the difference in distance between a-a/b-b pairs and 
a-b pairs. This distortion is quantified by the angle a 
between any three ions a, c, b; for a regular octahedron, 
a is a right angle, and so the distances a-a, b-b, a-b 
arc identical. In general, the split configuration can be 




(d) v 2 + (relaxed); A = 0.75/v^, / = 2 (e) v° (relaxed); A = 1.5/ y/v, f = 2 (f) v 1 " (relaxed); A = 1.75/y/v, / = 2 




(g) Split vacancy lattice schematic (h) v 2 (relaxed); A = 2.25/y/v, f = 2 



FIG. 6. (Color online) Contour-surface plots of the MLWFs most strongly associated with the defect centre for different charge 
states of the vacancy, calculated using the 256-atom BCC supercell and a norm-conserving pseudopotential; for computational 
efficiency we only include the T-point wavefunction in the wannierization procedure, since we expect a negligible difference in 
the qualitative features of the resulting Wannier functions. Each separate closed surface is an individual Wannier function. 
The amplitude A of the contour shown in each figure is given in the caption in terms of the primitive cell volume v. In general 
the Wannier functions also have small negative components on the backbonds (not pictured). / is the electronic occupancy of 
a single orbital. The wannierization procedure is performed using the WANNlEFtSO^i code (version 1.2). 

specified with only two parameters: the distortion angle six defect Wannier functions, each one a bond between 
a and the distance a-c between the central ion and any the central ion and one of its neighbours. This suggests 
of the ions forming the octahedral cage. that the c ion is forming six sp 3 d 2 orbitals which then 

bond to the dangling sp 3 orbitals of the a, b ions. The 

The wannierization of the occupied manifold produces 
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inwards relaxation of the ionic positions also shortens the 
bonds and reduces the distortion of the octahedron. This 
configuration is favoured since these bonds are shorter 
than the ones obtained by dimerization in the C2 V ar- 



rangement predicted by Watkins' model, and so arc closer 
in length to the bulk silicon bond. However, only in the 
case of the doubly negative vacancy are there enough 
electrons to fully occupy the six orbitals. 
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